Коммутирующие векторные поля

В.Э. Адлер. Классические интегрируемые системы, весенний семестр 2020

Если два векторных поля коммутируют, то отвечающие им динамические системы имеют совместное решение. В этой программе мы визуализируем это утверждение на численных примерах.

Немного теории

Пусть даны две конечномерные динамические системы

(1)	(Subscript[u, x])=a(u), a=((Subscript[a, 1]),...,(Subscript[a, n])),
(2)	(Subscript[u, y])=b(u),  b=((Subscript[b, 1]),...,(Subscript[b, n])),

где u=((Subscript[u, 1]),...,(Subscript[u, n])).  Можно ли построить функцию u(x,y), удовлетворяющую обеим системам? Легко понять, что, вообще говоря, это невозможно. Простейший контрпример доставляет задача восстановления скалярной функции v(x,y) по её частным производным:

(Subscript[v, x])=f(x,y),   (Subscript[v, y])=g(x,y)

(чтобы привести эту задачу к виду (1), (2), нужно ввести векторы u=(x,y,v), a=(1,0,f) и b=(0,1,g)). Ясно, что для того, чтобы такая функция v существовала (в некоторой односвязной области (x,y), в которой обе функции fи g определены и дифференцируемы), необходимо и достаточно, чтобы выполнялось равенство перекрёстных производных (Subscript[((Subscript[v, x])), y])=(Subscript[((Subscript[v, y])), x]), то есть

(Subscript[f, y])(x,y)=(Subscript[g, x])(x,y).

Если это выполнено, то ответ существует. Если нет – не существует. Всё просто.

Для систем общего вида тоже всё просто. В (1) правая часть a(u) – это векторное поле, заданное в фазовом пространстве переменных u или в некоторой его области U. Тогда для любой (гладкой) функции f(u) определена производная по x в силу системы (1), как результат применения линейного дифференциального оператора первого порядка

(Subscript[L, a])=(Subscript[a, 1])∂/∂(Subscript[u, 1])+...+(Subscript[a, n])∂/∂(Subscript[u, n]),

то есть

(Subscript[D, x])(f)=(Subscript[L, a])(f)=(Subscript[a, 1])∂f/∂(Subscript[u, 1])+...+(Subscript[a, n])∂f/∂(Subscript[u, n]).

В частности, компоненты самого векторного поля совпадают с производными от координатных функций, (Subscript[a, i])=(Subscript[L, a])((Subscript[u, i])). Оператор (Subscript[L, a]) можно отождествлять с самим векторным полем a. Часто, когда говорят о векторных полях, имеют в виду именно отвечающие им дифференцирования.

Критерий существования общего совместного решения систем (1), (2) заключается в перестановочности соответствующих операторов. Определим коммутатор двух векторных полей как оператор, действующий на произвольную функцию f по правилу

[(Subscript[L, a]),(Subscript[L, b])](f)=(Subscript[L, a])(Subscript[L, b])(f)-(Subscript[L, b])(Subscript[L, a])(f).

Отметим, что произведение (Subscript[L, a])(Subscript[L, b]) это дифференциальный оператор второго порядка, но в коммутаторе все вторые производные сокращаются, так что он снова является векторным полем (то есть, векторные поля образуют алгебру относительно операции коммутатора. Это важный пример алгебры Ли). Явная формула для коммутатора имеет вид

(3)	[(Subscript[L, a]),(Subscript[L, b])]=((Subscript[L, a])((Subscript[b, 1]))-(Subscript[L, b])((Subscript[a, 1])))∂/∂(Subscript[u, 1])+...+((Subscript[L, a])((Subscript[b, n]))-(Subscript[L, b])((Subscript[a, n])))∂/∂(Subscript[u, n]).

Если u(x,y) есть некоторое частное совместное решение систем (1), (2), то условие совпадения перекрёстных производных (Subscript[((Subscript[u, x])), y])=(Subscript[((Subscript[u, y])), x]) означает, что должно быть [(Subscript[L, a]),(Subscript[L, b])]((Subscript[u, i]))=0 для любой компоненты u. То есть, на рассматриваемом решении (когда мы заменяем везде u на заданную функцию u(x,y)) должно выполняться равенство

(4)	[(Subscript[L, a]),(Subscript[L, b])]=0.

Если же мы хотим, чтобы системы (1), (2) имели общее совместное решение, при каких угодно начальных условиях, то равенство (4) должно выполняться тождественно по u.

Перейдем к экспериментам.

Пара линейных систем с [A,B]=0

Рассмотрим такую пару:

(Subscript[u, x])=a=A u,   (Subscript[u, y])=b=B u,

где A и B постоянные матрицы. В этом случае

[(Subscript[D, x]),(Subscript[D, y])](u)=(Subscript[D, x])(B u)-(Subscript[D, y])(A u)=[B,A]u.

Таким образом, в данном случае коммутативность векторных полей это то же самое, что коммутативность матриц. Линейные уравнения можно решать аналитически, но мы не будем этого делать, ограничимся численной проверкой. Желающие могут вывести формулу для общего совместного решения в следующем примере.

Определим функцию

nsol[A,u0,{x0,x1}]  – численное решение задачи Коши (Subscript[u, x])=A u, u((Subscript[x, 0]))=(Subscript[u, 0]) при x∈[(Subscript[x, 0]),(Subscript[x, 1])].

nsol[{{a_,b_},{c_,d_}},{u10_,u20_},{x0_,x1_}]:=
NDSolve[
{u1'[x]==a u1[x]+b u2[x],
u2'[x]==c u1[x]+d u2[x],
u1[x0]==u10,
u2[x0]==u20},
{u1,u2},
{x,x0,x1}][[1]]

Зададим конкретные коммутирующие матрицы и начальное условие.

A={{1,3},
{-3,1}};

B={{1,0},
{0,1}};

A.B-B.A   (* коммутатор матриц *)

u0={1,1};

x0=y0=0;
x1=y1=2;
δ=ϵ=0.2;

Решим задачу Коши для первой системы от x=(Subscript[x, 0]) до (Subscript[x, 1]). 
Значения u в промежуточных точках, с некоторым интервалом δ, примем в качестве начальных условий для второй системы. 
Решим соответствующие задачи Коши для второй системы,  от y=(Subscript[y, 0]) до (Subscript[y, 1]). 
Построим фазовые траектории решений и на каждой из них отметим точки с интервалом ϵ.

(Вместо фазовых траекторий можно с тем же успехом рисовать интегральные кривые. Правда, так как у нас вектор имеет размерность 2 и есть еще 2 переменные x,y, то получатся кривые в 4-мерном пространстве. Поэтому, удобнее изображать решение в виде пары графиков для двух функций  (Subscript[u, 1])(x,y), (Subscript[u, 2])(x,y). Как это делать показано в разделе 4).

solA=nsol[A,u0,{x0,x1}];

solBA=Table[nsol[B,{u1[x0+i δ],u2[x0+i δ]}/.solA,{y0,y1}],{i,0,(x1-x0)/δ}];

pointsBA=Table[{u1[y0+j ϵ],u2[y0+j ϵ]}/.solBA[[i+1]],
{i,0,(x1-x0)/δ},{j,0,(y1-y0)/ϵ}];

grBA=Show[
{ParametricPlot[{u1[x]/.solA,u2[x]/.solA},{x,x0,x1},PlotStyle->Cyan],
ParametricPlot[
Table[{u1[y]/.solBA[[i+1]],u2[y]/.solBA[[i+1]]},{i,0,(x1-x0)/δ}],
{y,y0,y1},PlotStyle->Blue],
Graphics[{Red,PointSize[0.015],Point[#]&/@pointsBA}]
},
PlotRange->10{{-1,1},{-1,1}},
AspectRatio->Automatic,
AxesOrigin->{0,0},
AxesLabel->{"(Subscript[u, 1])","(Subscript[u, 2])"},
ImageSize->300
]

Сделаем все то же самое в другом порядке.

solB=nsol[B,u0,{x0,x1}];

solAB=Table[nsol[A,{u1[y0+j ϵ],u2[y0+j ϵ]}/.solB,{x0,x1}],{j,0,(y1-y0)/ϵ}];

pointsAB=Table[{u1[x0+i δ],u2[x0+i δ]}/.solAB[[j+1]],
{i,0,(x1-x0)/δ},{j,0,(y1-y0)/ϵ}];

grAB=Show[
{ParametricPlot[{u1[y]/.solB,u2[y]/.solB},{y,y0,y1},PlotStyle->Blue],
ParametricPlot[
Table[{u1[x]/.solAB[[j+1]],u2[x]/.solAB[[j+1]]},{j,0,(y1-y0)/ϵ}],
{x,x0,x1},PlotStyle->Cyan],
Graphics[{Red,PointSize[0.015],Point[#]&/@pointsAB}]
},
PlotRange->10{{-1,1},{-1,1}},
AspectRatio->Automatic,
AxesOrigin->{0,0},
AxesLabel->{"(Subscript[u, 1])","(Subscript[u, 2])"},
ImageSize->300
]

Убедимся, что отмеченные точки совпадают, в пределах точности вычислений, и совместим графики.

Max[Abs[Flatten[pointsBA-pointsAB]]]
Show[grAB,grBA]

Точки, вычисленные двумя разными способами, совпадают. Можно сначала сдвинуться на δ вдоль фазовой траектории первой системы и затем на ϵ вдоль фазовой траектории второй. Результат будет тот же, если эти сдвиги переставить. Это и означает совместность двух систем.

Пара линейных систем с [A,B]≠0

Скопируем предыдущий раздел и изменим матрицы так, чтобы они не коммутировали.

A={{1,3},
{-3,1}};

B={{1,0},
{0,2}};

A.B-B.A  (* коммутатор матриц *)

u0={1,1};

x0=y0=0;
x1=y1=2;
δ=ϵ=0.2;

solA=nsol[A,u0,{x0,x1}];

solBA=Table[nsol[B,{u1[x0+i δ],u2[x0+i δ]}/.solA,{y0,y1}],{i,0,(x1-x0)/δ}];

pointsBA=Table[{u1[y0+j ϵ],u2[y0+j ϵ]}/.solBA[[i+1]],
{i,0,(x1-x0)/δ},{j,0,(y1-y0)/ϵ}];

grBA=Show[
{ParametricPlot[{u1[x]/.solA,u2[x]/.solA},{x,x0,x1},PlotStyle->Cyan],
ParametricPlot[
Table[{u1[y]/.solBA[[i+1]],u2[y]/.solBA[[i+1]]},{i,0,(x1-x0)/δ}],
{y,y0,y1},PlotStyle->Blue],
Graphics[{Red,PointSize[0.015],Point[#]&/@pointsBA}]
},
PlotRange->10{{-1,1},{-1,1}},
AspectRatio->Automatic,
AxesOrigin->{0,0},
AxesLabel->{"(Subscript[u, 1])","(Subscript[u, 2])"},
ImageSize->300
]

solB=nsol[B,u0,{x0,x1}];

solAB=Table[nsol[A,{u1[y0+j ϵ],u2[y0+j ϵ]}/.solB,{x0,x1}],{j,0,(y1-y0)/ϵ}];

pointsAB=Table[{u1[x0+i δ],u2[x0+i δ]}/.solAB[[j+1]],
{i,0,(x1-x0)/δ},{j,0,(y1-y0)/ϵ}];

grAB=Show[
{ParametricPlot[{u1[y]/.solB,u2[y]/.solB},{y,y0,y1},PlotStyle->Blue],
ParametricPlot[
Table[{u1[x]/.solAB[[j+1]],u2[x]/.solAB[[j+1]]},{j,0,(y1-y0)/ϵ}],
{x,x0,x1},PlotStyle->Cyan],
Graphics[{Red,PointSize[0.015],Point[#]&/@pointsAB}]
},
PlotRange->10{{-1,1},{-1,1}},
AspectRatio->Automatic,
AxesOrigin->{0,0},
AxesLabel->{"(Subscript[u, 1])","(Subscript[u, 2])"},
ImageSize->300
]

Max[Abs[Flatten[pointsBA-pointsAB]]]
Show[grAB,grBA]

Сетка расползлась, красные точки не попадают в узлы. Теперь не всё равно, в каком порядке ходить по траекториям. Системы не совместны.

Нелинейный пример

Здесь мы сделаем то же самое для нелинейного уравнения, и заодно проиллюстрируем ещё два момента: 
	– как работать в том случае, если вместо динамической системы (1) задано ОДУ порядка n;
	– как выводить график решения вместо фазовых траекторий.

Задача. Пусть обозначено (Subscript[u, n])=(∂_x^n)(u). Проверить, символически и численно, что ОДУ

(5)	(Subscript[u, 4])+10u (Subscript[u, 2])+5(Subscript[u, 1]^2)+10(u^3)+(Subscript[c, 1]) ((Subscript[u, 2])+3(u^2))+(Subscript[c, 2])u+(Subscript[c, 3])=0,

где (Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]) произвольные постоянные, совместно с уравнением КдФ

(6)  	(Subscript[u, t])=(Subscript[u, 3])+6 u (Subscript[u, 1]).

Если это действительно так, построить график какого-нибудь решения u(x,t) удовлетворяющего обоим уравнениям.

Решение. Перепишем оба уравнения в виде динамических систем. Так как 4-я производная выражается из ОДУ (5) через младшие, то динамическими переменными являются (Subscript[u, 0]),(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3]) (для единообразия будем писать (Subscript[u, 0]) вместо u). Сразу можно написать векторное поле (Subscript[D, x]) (точнее – то, как оно действует на произвольную функцию от динамических переменных):

dx[f_]:=u[1]D[f,u[0]]+
u[2]D[f,u[1]]+
u[3]D[f,u[2]]-
(10u[0]u[2]+5((u[1])^2)+10((u[0])^3)+c1(u[2]+3((u[0])^2))+c2 u[0]+c3)D[f,u[3]]

С ОДУ справились. С КдФ немного посложнее. Это уравнение задаёт правило дифференцирования по t только для переменной (Subscript[u, 0]), а нам нужно ещё для (Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3]). Значит, нужно продолжить дифференцирование (Subscript[D, t]) на эти производные, пользуясь тем, что (Subscript[D, x]) уже определено и полагая, по определению, что (Subscript[D, t])((Subscript[u, k]))=(Subscript[D, x]^k)((Subscript[u, t])). Имеем:

ut=u[3]+6u[0]u[1];

dt[f_]:=ut D[f,u[0]]+
dx[ut]D[f,u[1]]+
dx[dx[ut]]D[f,u[2]]+
dx[dx[dx[ut]]]D[f,u[3]]

Можно посмотреть, что за векторное поле получилось. Вот его компоненты:

dt[u[0]]
Expand[dt[u[1]]]
Expand[dt[u[2]]]
Expand[dt[u[3]]]

Итак, построили два векторных поля. Считаем коммутатор. Для этого нужно вычислить его значения на координатных функциях (= динамических переменных), то есть, [(Subscript[D, x]),(Subscript[D, t])]((Subscript[u, k])), k=0,1,2,3.

Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

Всё сошлось. На самом деле, ясно, что равенство нулю коммутатора на первых трёх компонентах выполняется автоматически – в силу построения дифференцирования (Subscript[D, t]) по формуле (Subscript[D, t])((Subscript[u, k]))=(Subscript[D, x]^k)((Subscript[u, t])). Для примера, испортим что-то в уравнениях. Тогда коммутативность нарушится, но только на последней компоненте. Так что, можно было бы проверять только её, но лишний контроль не мешает.

(* портим ut *)
ut=u[3]+66u[0]u[1];
Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

(* откатываем назад *)
ut=u[3]+6u[0]u[1];
Table[Expand[dx[dt[u[k]]]-dt[dx[u[k]]]],{k,0,3}]

Переходим к численному счёту. Определяем решение задачи Коши по x.

nsolx[{u00_,u10_,u20_,u30_},{c1_,c2_,c3_},{x0_,x1_}]:=
NDSolve[
{u[0]'[x]==u[1][x],
u[1]'[x]==u[2][x],
u[2]'[x]==u[3][x],
u[3]'[x]==-(10u[0][x]u[2][x]+5((u[1][x])^2)+10((u[0][x])^3)+c1(u[2][x]+3((u[0][x])^2))+c2 u[0][x]+c3),
u[0][x0]==u00,
u[1][x0]==u10,
u[2][x0]==u20,
u[3][x0]==u30},
{u[0],u[1],u[2],u[3]},
{x,x0,x1}][[1]]

Определяем решение задачи Коши по t. Нам нужно добавить аргумент t к тем обозначениям, что у нас были. Сделаем это в следующей ячейке, а потом скопируем (хотя такая практика не есть хорошо. Для более серьезных задач нужен какой-то другой способ).

dt[u[0]]/.u[k_]:>u[k][t]
Expand[dt[u[1]]]/.u[k_]:>u[k][t]
Expand[dt[u[2]]]/.u[k_]:>u[k][t]
Expand[dt[u[3]]]/.u[k_]:>u[k][t]

nsolt[{u00_,u10_,u20_,u30_},{c1_,c2_,c3_},{t0_,t1_}]:=
NDSolve[
{u[0]'[t]==6 u[0][t] u[1][t]+u[3][t],
u[1]'[t]==-c3-c2 u[0][t]-3 c1 ((u[0][t])^2)-10 ((u[0][t])^3)+((u[1][t])^2)-c1 u[2][t]-4 u[0][t] u[2][t],
u[2]'[t]==-c2 u[1][t]-6 c1 u[0][t] u[1][t]-30 ((u[0][t])^2) u[1][t]-2 u[1][t] u[2][t]-c1 u[3][t]-4 u[0][t] u[3][t],
u[3]'[t]==c1 c3+c1 c2 u[0][t]+4 c3 u[0][t]+3 (c1^2) ((u[0][t])^2)+4 c2 ((u[0][t])^2)+22 c1 ((u[0][t])^3)+40 ((u[0][t])^4)-c1 ((u[1][t])^2)-40 u[0][t] ((u[1][t])^2)+(c1^2) u[2][t]-c2 u[2][t]+8 c1 u[0][t] u[2][t]+10 ((u[0][t])^2) u[2][t]-2 ((u[2][t])^2)-6 u[1][t] u[3][t],
u[0][t0]==u00,
u[1][t0]==u10,
u[2][t0]==u20,
u[3][t0]==u30},
{u[0],u[1],u[2],u[3]},
{t,t0,t1}][[1]]

Дальше – всё почти как в предыдущих примерах. Решаем задачу Коши для (5) от x=(Subscript[x, 0]) до (Subscript[x, 1]). Значения вектора ((Subscript[u, 0]),(Subscript[u, 1]),(Subscript[u, 2]),(Subscript[u, 3])) в промежуточных точках, с некоторым шагом δ, принимаем в качестве начальных условий для (6). Решаем соответствующие задачи Коши от t=(Subscript[t, 0]) до (Subscript[t, 1]). Потом делаем то же самое в другом порядке, с шагом ϵ по t. Сравниваем получившиеся сетки. Если значения совпадают – значит функция u(x,t) определена корректно, строим её 3D график по полученной сетке.

Следует отметить, что в этой задаче многие решения имеют полюсные особенности. Ниже, значения параметров и начальных условий подобраны так, чтобы их не было.

x0=0;
x1=55;
Mx=200; (* число точек в сетке по x *)

t0=0;
t1=25;
Mt=200; (* число точек в сетке по t *)

cc={0.15,-0.7,0}; (* параметры *)

U0={0.1,-0.1,-0.1,0.125};  (* начальные данные при x0, t0 *)

(* решаем систему по t; решаем систему по x для последовательности 
значений t[0]= t0, ..., t[Mt]=t1 и сохраняем результат в виде 2D массива *)

nst=nsolt[U0,cc,{t0,t1}];

nstx2=Table[
ns=nsolx[{u[0][t],u[1][t],u[2][t],u[3][t]}/.t->t0+j(t1-t0)/Mt/.nst,cc,{x0,x1}];
Table[u[0][x0+i(x1-x0)/Mx]/.ns,{i,0,Mx}],{j,0,Mt}];

(* делаем то же самое, меняя ролями x и t *)
nsx=nsolx[U0,cc,{x0,x1}];

nsxt2=Table[
ns=nsolt[{u[0][x],u[1][x],u[2][x],u[3][x]}/.x->x0+i(x1-x0)/Mx/.nsx,cc,{t0,t1}];
Table[u[0][t0+j(t1-t0)/Mt]/.ns,{j,0,Mt}],{i,0,Mx}];

(* сравниваем две сетки *)
Max[Abs[nstx2-Transpose[nsxt2]]]

Значения, вычисленные двумя способами, совпали, что и подтверждает коммутативность. Теперь можно использовать полученный массив точек, чтобы построить график решения.

ListPlot3D[nstx2,
DataRange->{{x0,x1},{t0,t1}},
PlotRange->{-1,2},
Mesh->None,
AxesLabel->Evaluate[Style[#,Italic,24]&/@{"x","t","u"}],
ImageSize->600,
BaseStyle->{FontFamily->"Times",FontSize->16}
]

Пояснение. То, что у нас получилось это пример так называемого двухзонного, или двухфазного решения КдФ. Уравнение (5) это один раз проинтегрированное стационарное уравнение для высшей симметрии КдФ пятого порядка плюс младшие симметрии с константами (Subscript[c, i]). Мы проверили, на примере, что такие стационарные уравнения приводят к паре коммутирующих векторных полей, что и позволяет использовать их для построения решений.